require(plm)
require(runjags)
require(coda)
require(wfe)

#### data preparation ####
JSP.data.Study1 <- read.csv("JSP_data_Study1.csv")

# exclude periods with no speech
# PS = plenary session, BC = Budget Committee, OC = other committees
JSP.personnel.period.date.Study1 <- as.Date(c("1979-01-20", "1980-02-10", "1980-12-03", 
                                              "1982-02-06", "1982-12-17", "1983-09-07", 
                                              "1984-02-28", "1985-01-19", "1985-12-18", 
                                              "1986-01-22", "1986-09-08", "1987-01-24", 
                                              "1988-02-13", "1989-01-25", "1990-04-05", 
                                              "1991-02-01", "1991-07-31", "1991-12-21", 
                                              "1993-01-09", "1993-09-25", "1994-01-12", 
                                              "1994-06-30", "1994-09-03", "1995-05-27", 
                                              "1995-08-08", "1995-09-21", "1996-01-19", 
                                              "1996-10-20"))

JSP.PS.data.Study1 <- JSP.BC.data.Study1 <- JSP.OC.data.Study1 <- JSP.data.Study1
for (i in 1:(length(JSP.personnel.period.date.Study1) - 1)) {
  temp.data <- subset(JSP.data.Study1, time == i)
  if (sum(temp.data$PS.char > 0) == 0) {
    JSP.PS.data.Study1 <- subset(JSP.PS.data.Study1, time != i)
  }
  if (sum(temp.data$BC.char > 0) == 0) {
    JSP.BC.data.Study1 <- subset(JSP.BC.data.Study1, time != i)
  }
  if (sum(temp.data$OC.char > 0) == 0) {
    JSP.OC.data.Study1 <- subset(JSP.OC.data.Study1, time != i)
  }
}

# convert data frames to the pdata.frame class of the plm package
JSP.PS.pdata.Study1 <- pdata.frame(JSP.PS.data.Study1, c("name_jp", "time"))
JSP.BC.pdata.Study1 <- pdata.frame(JSP.BC.data.Study1, c("name_jp", "time"))
JSP.OC.pdata.Study1 <- pdata.frame(JSP.OC.data.Study1, c("name_jp", "time"))

#### main analyses (Online Appendix E.2) ####
## SNTV period
Study1.JSP.PS <- plm(PS.char ~ leader + chair + minister + 
                       age + log.term + vote.share + 
                       E.37 + E.38 + E.39 + E.40, 
                     data = JSP.PS.pdata.Study1, 
                     model = "within")
round(summary(Study1.JSP.PS, 
              vcov = vcovHC(Study1.JSP.PS, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)
nrow(JSP.PS.pdata.Study1)  # number of observations
length(unique(JSP.PS.pdata.Study1$name_jp))  # number of MPs
length(unique(JSP.PS.pdata.Study1$time))  # number of periods

#### analysis of committees (Online Appendix E.3) ####
## Budget Committee
Study1.JSP.BC <- plm(BC.char ~ leader + chair + minister + 
                       age + log.term + vote.share + 
                       E.37 + E.38 + E.39 + E.40, 
                     data = JSP.BC.pdata.Study1, 
                     model = "within")
round(summary(Study1.JSP.BC, 
              vcov = vcovHC(Study1.JSP.BC, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)
nrow(JSP.BC.pdata.Study1)  # number of observations
length(unique(JSP.BC.pdata.Study1$name_jp))  # number of MPs
length(unique(JSP.BC.pdata.Study1$time))  # number of periods

## other committees
Study1.JSP.OC <- plm(OC.char ~ leader + chair + minister + 
                       age + log.term + vote.share + 
                       E.37 + E.38 + E.39 + E.40, 
                     data = JSP.OC.pdata.Study1, 
                     model = "within")
round(summary(Study1.JSP.OC, 
              vcov = vcovHC(Study1.JSP.OC, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)
nrow(JSP.OC.pdata.Study1)  # number of observations
length(unique(JSP.OC.pdata.Study1$name_jp))  # number of MPs
length(unique(JSP.OC.pdata.Study1$time))  # number of periods

round(exp(Study1.JSP.OC$coefficients[1]), 3)  # substantive interpretation

#### robustness check: other dependent variables (Online Appendix E.4) ####
## number of speeches
Study1.JSP.PS.count <- plm(PS.count ~ leader + chair + minister + 
                             age + log.term + vote.share + 
                             E.37 + E.38 + E.39 + E.40, 
                           data = JSP.PS.pdata.Study1, 
                           model = "within")
round(summary(Study1.JSP.PS.count, 
              vcov = vcovHC(Study1.JSP.PS.count, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.JSP.BC.count <- plm(BC.count ~ leader + chair + minister + 
                             age + log.term + vote.share + 
                             E.37 + E.38 + E.39 + E.40, 
                           data = JSP.BC.pdata.Study1, 
                           model = "within")
round(summary(Study1.JSP.BC.count, 
              vcov = vcovHC(Study1.JSP.BC.count, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.JSP.OC.count <- plm(OC.count ~ leader + chair + minister + 
                             age + log.term + vote.share + 
                             E.37 + E.38 + E.39 + E.40, 
                           data = JSP.OC.pdata.Study1, 
                           model = "within")
round(summary(Study1.JSP.OC.count, 
              vcov = vcovHC(Study1.JSP.OC.count, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

## making one or more speeches
Study1.JSP.PS.bin <- plm(PS.bin ~ leader + chair + minister + 
                           age + log.term + vote.share + 
                           E.37 + E.38 + E.39 + E.40, 
                         data = JSP.PS.pdata.Study1, 
                         model = "within")
round(summary(Study1.JSP.PS.bin, 
              vcov = vcovHC(Study1.JSP.PS.bin, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.JSP.BC.bin <- plm(BC.bin ~ leader + chair + minister + 
                           age + log.term + vote.share + 
                           E.37 + E.38 + E.39 + E.40, 
                         data = JSP.BC.pdata.Study1, 
                         model = "within")
round(summary(Study1.JSP.BC.bin, 
              vcov = vcovHC(Study1.JSP.BC.bin, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.JSP.OC.bin <- plm(OC.bin ~ leader + chair + minister + 
                           age + log.term + vote.share + 
                           E.37 + E.38 + E.39 + E.40, 
                         data = JSP.OC.pdata.Study1, 
                         model = "within")
round(summary(Study1.JSP.OC.bin, 
              vcov = vcovHC(Study1.JSP.OC.bin, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

#### robustness check: censored regression (Online Appendix E.4) ####
# function to draw initial values
initial.fun.Tobit <- function(seed1, seed2, n.var, n.id) {
  set.seed(seed1)
  alpha <- rnorm(1)
  beta <- rnorm(n.var)
  gamma <- rnorm(n.id)
  sigma <- runif(1, 1, 2)
  tau <- runif(1, 1, 2)
  list(alpha = alpha, beta = beta, 
       gamma = gamma, sigma = sigma, tau = tau, 
       .RNG.name = "base::Wichmann-Hill", .RNG.seed = seed2)
}
# draw initial values
Study1.JSP.PS.Tobit.inits <- list()
for (i in 1:3) {
  Study1.JSP.PS.Tobit.inits[[i]] <- 
    initial.fun.Tobit(100 * i, 1000 * i, 12, 
                      max(as.numeric(factor(JSP.PS.data.Study1$name_jp))))
}
# MCMC sampling
Study1.JSP.PS.Tobit <- run.jags(model = "random_effects_Tobit_JAGS.R", 
                                monitor = c("beta", "alpha", "sigma", "tau"), 
                                data = list(y = ifelse(JSP.PS.data.Study1$PS.char > 0, 
                                                       JSP.PS.data.Study1$PS.char, NA), 
                                            z = (JSP.PS.data.Study1$PS.char > 0) * 1, 
                                            x = scale(as.matrix(JSP.PS.data.Study1[, c("leader", "chair", "minister", 
                                                                                       "female", "age", "dynastic", 
                                                                                       "log.term", "vote.share", 
                                                                                       "E.37", "E.38", "E.39", "E.40")]), 
                                                      scale = FALSE), 
                                            id = as.numeric(factor(JSP.PS.data.Study1$name_jp)), 
                                            N = nrow(JSP.PS.data.Study1), 
                                            n.var = 12, 
                                            n.id = max(as.numeric(factor(JSP.PS.data.Study1$name_jp)))), 
                                inits = Study1.JSP.PS.Tobit.inits, 
                                n.chains = 3, burnin = 5000, sample = 5000, modules = "glm", 
                                summarise = FALSE, thin = 1, method = "parallel")
# convergence check
print(which(gelman.diag(Study1.JSP.PS.Tobit, multivariate = FALSE)[[1]][, 1] > 1.1))
# posterior summary
round(summary(Study1.JSP.PS.Tobit), 3)
# posterior probability that the coefficient of party leader (beta[1]) is negative
round(mean(unlist(Study1.JSP.PS.Tobit$mcmc[, 1]) < 0), 3)

Study1.JSP.BC.Tobit.inits <- list()
for (i in 1:3) {
  Study1.JSP.BC.Tobit.inits[[i]] <- 
    initial.fun.Tobit(100 * i, 1000 * i, 12, 
                      max(as.numeric(factor(JSP.BC.data.Study1$name_jp))))
}
Study1.JSP.BC.Tobit <- run.jags(model = "random_effects_Tobit_JAGS.R", 
                                monitor = c("beta", "alpha", "sigma", "tau"), 
                                data = list(y = ifelse(JSP.BC.data.Study1$BC.char > 0, JSP.BC.data.Study1$BC.char, NA), 
                                            z = (JSP.BC.data.Study1$BC.char > 0) * 1, 
                                            x = scale(as.matrix(JSP.BC.data.Study1[, c("leader", "chair", "minister", 
                                                                                       "female", "age", "dynastic", 
                                                                                       "log.term", "vote.share", 
                                                                                       "E.37", "E.38", "E.39", "E.40")]), 
                                                      scale = FALSE), 
                                            id = as.numeric(factor(JSP.BC.data.Study1$name_jp)), 
                                            N = nrow(JSP.BC.data.Study1), 
                                            n.var = 12, 
                                            n.id = max(as.numeric(factor(JSP.BC.data.Study1$name_jp)))), 
                                inits = Study1.JSP.BC.Tobit.inits, 
                                n.chains = 3, burnin = 5000, sample = 5000, modules = "glm", 
                                summarise = FALSE, thin = 1, method = "parallel")
print(which(gelman.diag(Study1.JSP.BC.Tobit, multivariate = FALSE)[[1]][, 1] > 1.1))
round(summary(Study1.JSP.BC.Tobit), 3)
round(mean(unlist(Study1.JSP.BC.Tobit$mcmc[, 1]) < 0), 3)

Study1.JSP.OC.Tobit.inits <- list()
for (i in 1:3) {
  Study1.JSP.OC.Tobit.inits[[i]] <- 
    initial.fun.Tobit(100 * i, 1000 * i, 12, 
                      max(as.numeric(factor(JSP.OC.data.Study1$name_jp))))
}
Study1.JSP.OC.Tobit <- run.jags(model = "random_effects_Tobit_JAGS.R", 
                                monitor = c("beta", "alpha", "sigma", "tau"), 
                                data = list(y = ifelse(JSP.OC.data.Study1$OC.char > 0, JSP.OC.data.Study1$OC.char, NA), 
                                            z = (JSP.OC.data.Study1$OC.char > 0) * 1, 
                                            x = scale(as.matrix(JSP.OC.data.Study1[, c("leader", "chair", "minister", 
                                                                                       "female", "age", "dynastic", 
                                                                                       "log.term", "vote.share", 
                                                                                       "E.37", "E.38", "E.39", "E.40")]), 
                                                      scale = FALSE), 
                                            id = as.numeric(factor(JSP.OC.data.Study1$name_jp)), 
                                            N = nrow(JSP.OC.data.Study1), 
                                            n.var = 12, 
                                            n.id = max(as.numeric(factor(JSP.OC.data.Study1$name_jp)))), 
                                inits = Study1.JSP.OC.Tobit.inits, 
                                n.chains = 3, burnin = 5000, sample = 5000, modules = "glm", 
                                summarise = FALSE, thin = 1, method = "parallel")
print(which(gelman.diag(Study1.JSP.OC.Tobit, multivariate = FALSE)[[1]][, 1] > 1.1))
round(summary(Study1.JSP.OC.Tobit), 3)
round(mean(unlist(Study1.JSP.OC.Tobit$mcmc[, 1]) < 0), 3)

#### robustness check: weighted linear fixed-effect model (Online Appendix E.4) ####
Study1.JSP.PS.WFE <- wfe(PS.char ~ leader + chair + minister + 
                           age + log.term + vote.share + 
                           E.37 + E.38 + E.39 + E.40, 
                         data = JSP.PS.data.Study1, 
                         unit.index = "name_jp", time.index = "time", 
                         treat = "leader", estimator = "fd")
round(summary(Study1.JSP.PS.WFE)$coefficients, 3)

Study1.JSP.BC.WFE <- wfe(BC.char ~ leader + chair + minister + 
                           age + log.term + vote.share + 
                           E.37 + E.38 + E.39 + E.40, 
                         data = JSP.BC.data.Study1, 
                         unit.index = "name_jp", time.index = "time", 
                         treat = "leader", estimator = "fd")
round(summary(Study1.JSP.BC.WFE)$coefficients, 3)

Study1.JSP.OC.WFE <- wfe(OC.char ~ leader + chair + minister + 
                           age + log.term + vote.share + 
                           E.37 + E.38 + E.39 + E.40, 
                         data = JSP.OC.data.Study1, 
                         unit.index = "name_jp", time.index = "time", 
                         treat = "leader", estimator = "fd")
round(summary(Study1.JSP.OC.WFE)$coefficients, 3)